miχPODS: New baseline#

TODO

  • more intake usage

  • why does the stability diagram load multiple times

  • Is there an error in ONI phase labeling for the newer runs, presumably SST changed some.

  • Add diagnostic for convective “adjustment” or diff == 1

  • hvplot defaults

    • turn off scroll zoom

    • copy over presentation stuff

Notes#

This notebook compares:

  1. TAO

  2. “old” baseline with KD=1e-5, KV=2e-4 (NCAR/MOM6#238)

  3. old baseline + kpp.lmd.004 with KPP ν0=2.5, Ric=0.2, Ri0=0.5

  4. new baseline : KD=0, KV=0

  5. new baseline : kpp.lmd.004 with KPP ν0=2.5, Ric=0.2, Ri0=0.5

  6. new baseline : kpp.lmd.005 with KPP ν0=2.5, Ric=0.3, Ri0=0.5

Summary#

  1. Turning down the background visc by a factor of 40, (2e-4 → 5e-5 m2/s)

    1. sharpens the EUC relative to TAO

    2. changes the stability diagram for El-Nino warming phase.

  2. Using Ri_c=0.2, so shallower KPP surface layer, is key as Dan mentioned.

Catalog display#

Hide code cell source 
%load_ext rich

# MOM6 run catalog
catalog = {
    "baseline": (
        "baseline",
        "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods",
    ),
    "epbl": ("ePBL", "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.epbl.001.mixpods"),
    "kpp.lmd.002": (
        "KPP Ri0=0.5",
        "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.002.mixpods",
    ),
    "kpp.lmd.003": (
        "KPP Ri0=0.5, Ric=0.2,",
        "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.003.mixpods",
    ),
    "kpp.lmd.004": (
        "KPP ν0=2.5, Ric=0.2, Ri0=0.5",
        "gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods",
    ),
    "baseline.N150": (
        "baseline N=150",
        "gmom.e23.GJRAv3.TL319_t061_zstar_N150.baseline.mixpods",
    ),
    "kpp.lmd.004.N150": (
        "KPP ν0=2.5, Ric=0.2, Ri0=0.5, N=150",
        "gmom.e23.GJRAv3.TL319_t061_zstar_N150.kpp.lmd.004.mixpods",
    ),
    "new_baseline.hb": (
        "KD=0, KV=0",
        "gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.hb",
    ),
    "new_baseline.kpp.lmd.004": (
        "KPP ν0=2.5, Ric=0.2, Ri0=0.5",
        "gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.004.mixpods",
    ),
    "new_baseline.kpp.lmd.005": (
        "KPP ν0=2.5, Ri0=0.5",
        "gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.005.mixpods",
    ),
}
catalog

The rich extension is already loaded. To reload it, use:
  %reload_ext rich

{
    'baseline': ('baseline', 'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods'),
    'epbl': ('ePBL', 'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.epbl.001.mixpods'),
    'kpp.lmd.002': (
        'KPP Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.002.mixpods'
    ),
    'kpp.lmd.003': (
        'KPP Ri0=0.5, Ric=0.2,',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.003.mixpods'
    ),
    'kpp.lmd.004': (
        'KPP ν0=2.5, Ric=0.2, Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods'
    ),
    'baseline.N150': ('baseline N=150', 'gmom.e23.GJRAv3.TL319_t061_zstar_N150.baseline.mixpods'),
    'kpp.lmd.004.N150': (
        'KPP ν0=2.5, Ric=0.2, Ri0=0.5, N=150',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N150.kpp.lmd.004.mixpods'
    ),
    'new_baseline.hb': ('KD=0, KV=0', 'gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.hb'),
    'new_baseline.kpp.lmd.004': (
        'KPP ν0=2.5, Ric=0.2, Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.004.mixpods'
    ),
    'new_baseline.kpp.lmd.005': (
        'KPP ν0=2.5, Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.005.mixpods'
    )
}

catalog_sub = {
    k: catalog[k]
    for k in catalog.keys()
    if k == "kpp.lmd.004" or ("baseline" in k and "150" not in k)
}
display(catalog_sub)

{
    'baseline': ('baseline', 'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.001.mixpods'),
    'kpp.lmd.004': (
        'KPP ν0=2.5, Ric=0.2, Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.baseline.kpp.lmd.004.mixpods'
    ),
    'new_baseline.hb': ('KD=0, KV=0', 'gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.hb'),
    'new_baseline.kpp.lmd.004': (
        'KPP ν0=2.5, Ric=0.2, Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.004.mixpods'
    ),
    'new_baseline.kpp.lmd.005': (
        'KPP ν0=2.5, Ri0=0.5',
        'gmom.e23.GJRAv3.TL319_t061_zstar_N65.new_baseline.kpp.lmd.005.mixpods'
    )
}

Setup#

Hide code cell source 
%load_ext autoreload
%load_ext rich
%load_ext watermark

import math
import warnings

import cf_xarray as cfxr
import dcpy.datatree  # noqa
import distributed
import holoviews as hv
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import tqdm
from datatree import DataTree

import xarray as xr

%aimport pump
from pump import mixpods

hv.notebook_extension("bokeh")

cfxr.set_options(warn_on_missing_variables=False)
xr.set_options(keep_attrs=True, display_expand_data=False)

plt.style.use("bmh")
plt.rcParams["figure.dpi"] = 140


%watermark -iv
Hide code cell output 
The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
The rich extension is already loaded. To reload it, use:
  %reload_ext rich

matplotlib : 3.7.1
holoviews  : 1.15.4
numpy      : 1.23.5
distributed: 2023.3.0
pandas     : 1.5.3
pump       : 1.0+247.g1f1c5e1.dirty
cf_xarray  : 0.8.0
tqdm       : 4.65.0
xarray     : 2023.3.0
Hide code cell source 
if "client" in locals():
    client.close()
    del client
if "cluster" in locals():
    cluster.close()

import ncar_jobqueue

cluster = ncar_jobqueue.NCARCluster(
    local_directory="/local_scratch/pbs.$PBS_JOBID/dask/spill",
)
cluster.adapt(minimum_jobs=1, maximum_jobs=4)
client = distributed.Client(cluster)
client
Hide code cell output 
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/dask_jobqueue/core.py:20: FutureWarning: tmpfile is deprecated and will be removed in a future release. Please use dask.utils.tmpfile instead.
  from distributed.utils import tmpfile

Client

Client-019abd0a-d491-11ed-9153-3cecef1aca66

Connection method: Cluster object Cluster type: dask_jobqueue.PBSCluster
Dashboard: https://jupyterhub.hpc.ucar.edu/stable/user/dcherian/proxy/8787/status

Cluster Info

Read#

Subset Catalog, build tree#

%autoreload
datasets = {}
with warnings.catch_warnings():
    warnings.simplefilter("ignore", category=UserWarning)
    warnings.simplefilter("ignore", category=FutureWarning)
    for short_name, (long_name, folder) in tqdm.tqdm(catalog_sub.items()):
        datasets.update(
            {
                short_name: mixpods.load_mom6_sections(folder).assign_attrs(
                    {"title": long_name}
                )
            }
        )

100%|██████████| 5/5 [00:30<00:00,  6.10s/it]

tree = DataTree.from_dict(datasets)
tree

<xarray.DatasetView>
Dimensions:  ()
Data variables:
    *empty*

Modify#

tao_gridded = mixpods.load_tao()
tree["TAO"] = DataTree(tao_gridded)

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/dataset.py:255: UserWarning: The specified Dask chunks separate the stored chunks along dimension "depth" starting at index 42. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/dataset.py:255: UserWarning: The specified Dask chunks separate the stored chunks along dimension "time" starting at index 199726. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/dataset.py:255: UserWarning: The specified Dask chunks separate the stored chunks along dimension "longitude" starting at index 2. This could degrade performance. Instead, consider rechunking after loading.
  warnings.warn(

tree = tree.dc.reorder_nodes(["TAO", ...])

datasets["les"] = les["0.0.-140.oct.average.month"].to_dataset()

# Offset LES to work with slicing below
datasets["les"]["time"] = datasets["les"]["time"] + pd.Timedelta(days=25 * 365)

if "les" in tree:
    tree["les"] = tree["les"].isel(z=slice(-2))
    tree["les"]["KT"].attrs["standard_name"] = "ocean_vertical_heat_diffusivity"

if "micro" in locals():
    tree.update(micro)

Post-process catalog subset#

tree = tree.sel(time=slice("2000", "2017"))
tree

<xarray.DatasetView>
Dimensions:  ()
Data variables:
    *empty*
ref = tree["TAO"].ds.reset_coords(drop=True).cf.sel(Z=slice(-120, None))
counts = np.minimum(ref["S2"].cf.count("Z"), ref["N2T"].cf.count("Z")).load()


def calc_histograms(ds):
    ds = ds.copy()
    ds["tao_mask"] = counts.reindex(time=ds.time, method="nearest") > 5
    ds["tao_mask"].attrs = {
        "description": "True when there are more than 5 5-m T, u, v in TAO dataset"
    }
    # ds = ds.where(ds.tao_mask)
    return ds.update(mixpods.pdf_N2S2(ds))


tree = tree.map_over_subtree(calc_histograms)

daily composite#

%autoreload

dailies = tree.map_over_subtree(mixpods.daily_composites)
dailies

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out,

<xarray.DatasetView>
Dimensions:  ()
Data variables:
    *empty*

seasonal cycle#

tree = mixpods.persist_tree(tree)

Validate before continuing#

mixpods.validate_tree(tree)

Verify depth is normalized

for name, ds in datasets.items():
    (ds.cf["sea_water_x_velocity"].cf["Z"].plot(marker=".", ls="none", label=name))
plt.legend()

<matplotlib.legend.Legend>
_images/2372712bb196933da3d938a29d3f04e0e217fe29d04bb036dbceb08ab9ce12da.png

Mean profiles#

State Variables#

  1. A lot more shear, \(S, S^2\) just above the EUC when visc is turned down!

  2. We are very slightly lower on \(S^2\), \(N_T^2\) in the top 80m, compare TAO, kpp.lmd.004, new_baseline.kpp.lmd.004.

  3. Using the standard Ri_c=0.3, so deeper KPP surface layer, decreases \(S^2\), \(N_T^2\) in the top 60m. new_baseline.kpp.lmd.004 vs new_baseline.kpp.lmd.005

for context, Peters et al (1995) is interesting:

  1. Variability patterns at 50-350 m are distinctly different from the upper 50 m containing the diurnal cycle of mixing

  2. Large-scale shear of wavenumbers k < 0.01 cpm is associated with the slowly varying EUC and EIC. Large-scale is the dominant component of total shear above 50 m and in the thermostad, 170-270 m.

  3. Fine-scale shear, 0.01cpm < k < 0.5 cpm,provides the dominant component of total rms shear and exceeds the large-scale component in thick layersaroundthe coresof EUC and EIC, where Fr < 1, at 50-170 m and below 270 m.

  4. Variations of fine-scale shear cause variations in turbulent mixing; the large-scaleshear alone is a poor predictor of mixing.

  5. Lacking a model that links large-scale, fine-scale, and turbulent flow components,our service to equatorial modelers consists of describing general levels of eddy diffusivities and their variability patterns.

These models should not be resolving “finescale” shear, and the mixing is not well correlated with “large-scale shear”, so we need a “background” diffusivity/viscosity to make things work.

But see that std(\(S^2\)) is a lot stronger above the EUC core, relative to TAO for new_baseline.kpp.lmd.004 and new_baseline.kpp.lmd.005

%autoreload

S2 = mixpods.plot_profile_fill(tree, "S2", "S^2")
N2 = mixpods.plot_profile_fill(tree, "N2T", "N_T^2")
u = mixpods.plot_profile_fill(tree, "sea_water_x_velocity", "u")
T = mixpods.plot_profile_fill(tree, "sea_water_potential_temperature", "T")
Ri = mixpods.plot_median_Ri(tree)

(S2 + N2 + u + T).cols(2).opts(
    hv.opts.Curve(xlim=(-350, 0)),
    hv.opts.Layout(toolbar="left"),
    clone=True,
)

Ri

Turbulence Variables#

ENSO Phases#

Stability Diagram#

  1. Major change in La-Nina Warming for the new-baseline runs. We need to check ONI closely.

%autoreload

mixpods.plot_stability_diagram_by_dataset(tree, nrows=2)
_images/6840677fd142da5c4dd802da779c8f4282e6cf0500c6ea98ec4efc70cf8238b7.png

ε-Ri histograms#

mixpods.map_hvplot(
    lambda dt, name, muted: mixpods.plot_eps_ri_hist(
        dt["eps_ri"].load(), label=name, muted=muted
    ),
    tree.children,
).opts(show_legend=True).opts(hv.opts.Curve(ylim=(1e-8, 2e-6)))

Daily composites#

Hard to interpret! I think a lot of this is bias in KPP surface layer vs actively mixing layer in obs. We can write better diagnostics to check this (Moum et al, 2023)

  1. The 89m χpod comparison is quite interesting. Suggests we do need more background visc

%autoreload

mixpods.plot_daily_composites(dailies, ["eps", "KT"], logy=True, legend=False).opts(
    hv.opts.GridSpace(show_legend=True)
)

%autoreload

mixpods.plot_daily_composites(
    dailies, ["S2", "N2", "Rig_T"], logy=False, legend=False
).opts(hv.opts.GridSpace(show_legend=True, width=200), hv.opts.Overlay(frame_width=200))

Turbulence distributions#

handles = [
    mixpods.plot_distributions(tree, "chi", bins=np.linspace(-11, -4, 101), log=True),
    mixpods.plot_distributions(tree, "eps", bins=np.linspace(-11, -4, 101), log=True),
    mixpods.plot_distributions(
        tree, "ocean_vertical_heat_diffusivity", bins=np.linspace(-8, -1, 101), log=True
    ),
    # plot_distributions(tree, "Jq", bins=np.linspace(-1000, 0, 51), log=False),
    mixpods.plot_distributions(tree, "Rig_T", np.linspace(-0.5, 1.5, 61))
    * hv.VLine(0.25).opts(line_color="k"),
]

/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/computation.py:769: RuntimeWarning: invalid value encountered in log10
  result_data = func(*input_data)
/glade/u/home/dcherian/miniconda3/envs/pump/lib/python3.10/site-packages/xarray/core/computation.py:769: RuntimeWarning: invalid value encountered in log10
  result_data = func(*input_data)

hv.Layout(handles).opts("Overlay", frame_width=600).cols(2)

Compare boundary layer depth#

mixing_layer_depth_criteria = {
    "boundary_layer_depth": {"name": "KPPhbl|KPP_OBLdepth|ePBL_h_ML"},
}

hbl = (
    tree.drop_nodes("TAO")
    .dc.subset_nodes("KPP_OBLdepth")
    .dc.concatenate_nodes()
    .reset_coords(drop=True)
).load()

(
    # hbl.groupby("time.hour").mean().hvplot.line(by="node", flip_yaxis=True)
    hbl.groupby("time.hour").median().hvplot.line(by="node", flip_yaxis=True)
    + hbl.to_dataframe().hvplot.hist(
        by="node",
        bins=np.arange(0, 90, 1),
        normed=1,
        alpha=0.3,
        ylim=(0, 0.05),
        muted_alpha=0,
    )
).opts(hv.opts.Curve(invert_yaxis=True)).opts(mixpods.HV_TOOLS_OPTIONS).opts(
    mixpods.PRESENTATION_OPTS
)

Heat Budget#

f, ax = plt.subplots(
    2,
    math.ceil(len(tree) / 2),
    constrained_layout=True,
    squeeze=False,
    sharex=True,
    sharey=True,
    figsize=(10, 3),
)

for axx, (name, node) in zip(ax.flat, tree.children.items()):
    mixpods.plot_climo_heat_budget_1d(node.ds, mxldepth=-40, penetration="mom", ax=axx)
    axx.set_title(name)

dcpy.plots.clean_axes(ax)